clear; close all
qbit=1;showfig = 0;

Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
L = [0.04 0.2:0.2:3];

factor = 2.5;
x1 = 40e-1*factor;
y1 = 50e-1*ones(size(factor));
z1 = 90e-1*factor;
x2 = -30e-1*factor;
y2 = 50e-1*ones(size(factor));
z2 = 40e-1*factor;
TRX = zeros(6,length(factor));
TRX(1,:) = x1;
TRX(2,:) = y1;
TRX(3,:) = z1;
TRX(4,:) = x2;
TRX(5,:) = y2;
TRX(6,:) = z2;
load('TRX_std.mat','TRX','NTRX')
NTRX = 200;
dist = rms(TRX(1:3,:))*sqrt(3)+rms(TRX(4:6,:))*sqrt(3);

fc = 30e9;
c = 3e8;
BW_des = 1.5e9;
Ptx = 10;%dBm

cal_rate = 0;
if cal_rate

Rates = zeros(3,length(Ptx),NTRX);
Ratessub = Rates;
RatesGS = zeros(1,length(Ptx),NTRX);
for n = 1:length(L)
    len = L(n);
    parfor m = 1:NTRX
   Rates(:,n,m) = my_compute_rate(10,Nsub,len,TRX(:,m),fc,BW_des,Ptx,showfig); 
    Ratessub(1,n,m) = compute_rate_sub(10,Nsub,len,TRX(:,m),fc,BW_des,Ptx,showfig,[2 0]); 
    RatesGS(1,n,m) = compute_rate_GS(10,Nsub,len,TRX(:,m),fc,BW_des,Ptx,showfig); 
    end
    disp(n)
end
Rates = mean(Rates,3);
Ratessub = mean(Ratessub,3);
RatesGS = mean(RatesGS,3);
save('data/Rates_sum_rate_vs_D.mat','Rates',"Ratessub",'RatesGS','L')
else
    load('data/Rates_sum_rate_vs_D.mat');
end

figure(1)
plot(L,Rates(3,:),'--','LineWidth',1.5) % optimal
hold on
grid on
plot(L,Rates(2,:),'v-','LineWidth',1.5) % proposed
plot(L,RatesGS(1,:),'*-','LineWidth',1.5)
plot(L,Rates(1,:),'k^--','LineWidth',1.5); % classial 
plot(L,Ratessub(1,:),'+--','LineWidth',1.5); % classial 
xlabel('RIS side length D (m)')
ylabel('Average achievable rate (bit/s)')
% legend('Optimal beamforming','Proposed method','Proposed method, b = 3' ...
%     ,'Proposed method, b = 2','Proposed method, b = 1','Classical beamforming','Classical beamforming, b = 1' ...
%     ,'Location','southeast')
legend('Optimal beamforming','FZ-SPM','FZ-GSA','Classical beamforming','VSA-based beamforming',...
    'Location','southeast', 'FontSize', 14)
hold off
set(get(gca, 'Legend'), 'Box', 'off'); 
set(gca, 'FontSize', 16);
set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小
set(gca,'FontName','Times New Roman')
set(0, 'DefaultAxesFontName', 'Times New Roman');
saveas(gcf,'fig/Fig_sum_rate_vs_D.fig')
saveas(gcf,'fig/Fig_sum_rate_vs_D.png')

figure(2)
plot(L,Rates(3,:),'--','LineWidth',1.5) % optimal
hold on
% plot(L,Rates(2,:),'v-','LineWidth',1.5) % proposed
plot(L,Rates(1,:),'k^--','LineWidth',1.5); % classial 
xlabel('RIS side length D (m)')
ylabel('Average rate per subcarrier (bit/s)')
% legend('Optimal beamforming','Proposed method','Proposed method, b = 3' ...
%     ,'Proposed method, b = 2','Proposed method, b = 1','Classical beamforming','Classical beamforming, b = 1' ...
%     ,'Location','southeast')
legend('Optimal beamforming','Classical beamforming',...
    'Location','southeast');
hold off
